This post is to simply introduce how to make a map in R environment, or more officially, geographic/spatial visualisation. However, I just found numerous methods to do this in R, depending on which data type you have. So this post will not introduce the basic of R and GIS. After all, I am also not an expert in this.
To make a map, at least two types of data are needed: * Geographic Information System (GIS) data for position * Research data per se for value
Sometimes these two are combined to a unified data (e.g., population
data as an attribute in a shapefile). However, in most cases for those
who are not studying geography (like me), we have to join the research
data (in csv, netCDF) to GIS fields based on shared “key/id column”. The
relevant GIS data probably can be found in third-party packages such as
maps or mapdata, but for more demands like
custom resolution, one have to download in other sources.
There are two methods to plot a csv data frame on a map:
plot layers using ggplot2
transform to sf
Because the latter one is same to sections below, so I will introduce
ggplot2, which is very (but not limitedly) proficient at
analysing data frames.
library(ggplot2)
library(maps) # outline data for continents,countries
library(mapdata) # extra map datasets
library(mapproj) # map projection
#library(maptools)
# a quickest method to plot a global map
# map("world")
# Here I create a simple dataframe, which can be replaced to your data by read.csv()
toy_df <- data.frame(city = c("Bristol", "Shangrao"),
long = c(51.45, 28.45),
lat = c(-2.59, 117.94))
world <- map_data("world") # a ggplot2 function to get a world data.frame,see also borders
# notice x is lat and y is long, which differs normal impression
ggplot(world) +
geom_polygon(aes(x = long, y = lat, group=group)) +
geom_label(data= toy_df, aes(lat, long, label=city)) +
coord_map("mollweide", xlim=c(-180,180)) + # change the projection
theme_minimal()
Actuallyggfortify provides the autoplot
function, which is a quicker wrapper for ggplot2. So we can do the same
by using less command.
library(ggfortify)
world <- map('world', plot = FALSE, fill = TRUE)
autoplot(world, geom = 'polygon') +
geom_point(data = toy_df, aes(lat, long),
color="red", size=2)
There is another geom type: geom_map in ggplot. The difference between geom_map and geom_polygon is that the latter just plot the polygon (or position) and does not care about the value, but the geom_map contains both position and value information. Therefore, there must be a id columns shared by both polygon and value data frames, but no need to add lat/long in aes().
An example
fr_map <- map_data("france")
fr_data <- data.frame(region = unique(fr_map$region),
value = rnorm(length(unique(fr_map$region))))
ggplot(fr_data, aes(map_id = region)) +
geom_map(aes(fill = value), map=fr_map) +
expand_limits(x = fr_map$long, y = fr_map$lat) +
coord_fixed()
shp file can be read in R through sf or
rgal, and then plot using ggplot .
#from https://r-spatial.org/r/2018/10/25/ggplot2-sf.htm
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rnaturalearth) #another public map data package
#library(rnaturalearthdata)
it_data =st_as_sf(map("italy", plot = FALSE, fill = TRUE)) #polygons
it_sf <- ne_countries(scale = "medium", country="italy",returnclass = "sf") #an outline
ggplot(it_sf)+geom_sf()+geom_sf(data = it_data, fill = NA)
#and you can add more geom_sf() based on data in your hand
#plot(uk_sf$geometry) #uk_sf$geometry equals st_geomety(uk_sf)
tmap has a similar layer-based syntax to plot map, but
should be more professional and specilised than
ggplot2.
#from https://geocompr.robinlovelace.net/adv-map.html
library(spData)
#library(spDataLarge)
library(tmap)
library(sf)
world <- spData::world
toy_sf <- st_as_sf(toy_df, coords = c("lat", "long"), crs=4326) #WGS84
tm_shape(world)+tm_borders()+tm_fill()+
tm_shape(toy_sf) + tm_dots(size=.5)
ggmap is another friendly package for people from other
research background, it provides a function get_map() to
download some map data from google, open street map and so on.
library(ggmap)
#downloading raster data by get_map(). REQUIRE GOOGLE API
gmap <- get_stamenmap(bbox = c(left = -10, bottom = 50, right = 5, top = 60),
zoom = 6,
maptype='watercolor',
messaging = FALSE)
ggmap(gmap) + theme_void()
netCDF data can be transformed to raster format.
library(ncdf4)
library(terra)
library(rasterVis)
# an example file from UCAR
# https://www.unidata.ucar.edu/software/netcdf/examples/files.html
nc <- nc_open("tos_O1_2001-2002.nc")
#names(nc$var) # print all variables
sst_K <- ncvar_get(nc, varid="tos") #in Kelvin
sst <- sst_K - 273.15 #in Celsius
#dim(sst) # check dimension
#image(sst[,,1])
sst_raster <- flip(rast(t(sst[,,1])))
plot(sst_raster)
We can also use ggplot to do the plotting. But notice that the array from netCDF is not containing CRS data, that’s why we have x-y axis with 0-1 range.
library(ggplot2)
#to use geom_raster/geom_tile, we should transfer raster to a data.frame
sst_df <- sst_raster |> terra::as.data.frame(xy=TRUE) #pipe requires R > 4.1, otherwise use %>% in tidyverse
names(sst_df)[3] <- "sst"
ggplot(sst_df) + geom_raster(aes(x=x, y=y, fill=sst)) + theme_light() + scale_fill_viridis_c()
# what's the difference between |> and %>%?
#https://stackoverflow.com/questions/67633022/what-are-the-differences-between-rs-new-native-pipe-and-the-magrittr-pipe
As for the difference between geom_tile,
geom_raster, let’s see the offical answer:
geom_rect() and geom_tile() do the same thing, but are parameterised differently: geom_rect() uses the locations of the four corners (xmin, xmax, ymin and ymax), while geom_tile() uses the center of the tile and its size (x, y, width, height). geom_raster() is a high performance special case for when all the tiles are the same size.
rworldmap for extra intersting datacowplot for subplot added on ggplot objectsgganimate for animated plotleaflet for interactive map